function patient = Generate_4DCT_toolbox(patient, use_average_image, num_phases )
    
    if nargin < 2
        %         use_average_image = 0; %default - don't used the averaged image to deform (use the reference instead).
        prompt = 'Would you like to deform the averaged image [1] or the referece image [0] ?';
        use_average_image = input(prompt);
    end
    
    if nargin<3
        num_phases = 8; %Default number of phases = 8 (as in current clinical Siemens 4D-CT);
    end
%     patient.num_phases_4dct = num_phases;
    if use_average_image>0;
        load([patient.folder '/img_average'])
        deform_image = img_average;
    else
        deform_image = patient.static;
    end
    
    patient.folder_4DCT = [patient.model_folder '/Psuedo_4DCT'];
    mkdir(patient.folder_4DCT)
    cd(patient.folder_4DCT)
    dirinfo = dir('*phase*');
    dirinfo(~[dirinfo.isdir]) = [];  %remove non-directories
    
    str = 'Y'; %default is to continue
    if size(dirinfo,1)== num_phases;
        prompt = 'Psuedo 4DCT exists. \n                Do you want to re-generate 4DCT? Y [Hit Return for N]: ';
        str = input(prompt,'s');
        if isempty(str)
            str = 'N';
        end
    end
    if strcmpi('Y', str)>0;str = 'Y';end
     if strcmpi('N', str)>0;str = 'N';end
   
    if str == 'N'
        
    else
        
        
        step = 1;
        ii = 1:step:patient.dim(1);
        jj = 1:step:patient.dim(2);
        kk = 1:step:patient.dim(3);
        
        
        % img_average = static;
        % img_average(body == 0) = min(static(:));
        %% ------------------
        %Generate a 'representative' breath, from the bellows_voltage signal;
        
        [bellows_volt_rep, flow_rep, original_bellows]=gernerate_representative_breath_toolbox(patient);%folder_static, drift, pt)
        [bellows_orig, flow_orig]=calculate_volume_flow_fit_LongBreath_toolbox(original_bellows, 0.01);%CHECK RATE;
        figure;plot(bellows_orig, flow_orig,'b')
        hold on;
        plot(bellows_volt_rep, flow_rep,'r', 'LineWidth', 2)
        
        
        %double the breath, incase it doesn't cover the required phases;
        breath_rep = [bellows_volt_rep];% bellows_volt_rep];
        flow_rep = [flow_rep];% flow_rep];
        
        
%             v5 = -prctile(-patient.bellows_volt_drifted(:),5)
                        v5 = -prctile(-bellows_orig(:),5)

%             v5 = -prctile(-max(patient.bellows_volt_drifted),5)
        %
%             v95 = -prctile(-patient.bellows_volt_drifted(:),95)
                        v95 = -prctile(-bellows_orig(:),95)

%             v95 = -prctile(-min(patient.bellows_volt_drifted),95)
        %
%             f5 = prctile(patient.flow_drifted(:),5)
%             f85 = prctile(patient.flow_drifted(:),85)
        %
        %     figure;hold on
        %     plot(patient.bellows_volt_drifted,'b');
        %     plot([1 patient.dim(3)], [v5 v5],'g','LineWidth',2);
        %     plot([1 patient.dim(3)], [v95 v95],'r','LineWidth',2);
        %     hold off
             rep_breath_tmp = abs(v5-v95)*(bellows_volt_rep-min(bellows_volt_rep))./(max(bellows_volt_rep)-min(bellows_volt_rep))+v95;
%              rep_flow_tmp = abs(f5-f85)*(flow_rep-min(flow_rep))./(max(flow_rep)-min(flow_rep))-f85;
        %
        %     figure;
        %     plot(rep_breath_tmp)
        %     hold on
        %     plot([1 length(rep_breath_tmp)], [v5 v5],'g','LineWidth',2);
        %     plot([1 length(rep_breath_tmp)], [v85 v85],'r','LineWidth',2);
        
        breath_rep = rep_breath_tmp
        phases=[0:25:100]
        figure;plot(breath_rep,'.')
        volt_prctiles = -prctile(-breath_rep,phases)
        volt_prctiles = [volt_prctiles volt_prctiles([4,3,2])]
        phases = [phases];% phases([4,3,2])]
        
        hold on;
        plot([1 length(breath_rep)], [volt_prctiles' volt_prctiles'])
        for p = 1:length(phases);
            [~, point_x(p)] = min(abs(breath_rep(flow_rep>0)-volt_prctiles(p)));
        end
        for p = length(phases)+1:num_phases;
            [~, point_x(p)] = min(abs(breath_rep(flow_rep<=0)-volt_prctiles(p)));
        end
        point_x = sort(point_x);
        %point_x = indexs
        %     plot(indexs, breath_rep(indexs),'ro')
             title('Choose 8 phase points at the interception of the breath and the 8%iles;')
             [point_x,point_y] = ginput(8)
        
        for p = 1:length(phases);
            [~, indexs] = sort(abs(breath_rep-volt_prctiles(p)),'ascend');
            indexs = indexs(1:10)
            [~, pp] = min(abs(point_x(p)-indexs))
            point(p) = indexs(pp)
        end
        for p = length(phases)+1:num_phases;
            [~, indexs] = sort(abs(breath_rep-volt_prctiles(p)),'ascend');
            indexs = indexs(1:10)
            [~, pp] = min(abs(point_x(p)-indexs))
            point(p) = indexs(pp)
        end
        
        plot(point,breath_rep(point),'ro','LineWidth',2)
        %     close all
        figure(100);hold on
        plot(bellows_orig, flow_orig)
        plot(breath_rep, flow_rep,'r','LineWidth',4)
        hold on
        plot([volt_prctiles' volt_prctiles'], [min(flow_rep) max(flow_rep)],'LineWidth',4)        
        plot(breath_rep(point), flow_rep(point),'ro','LineWidth',4)

        xlabel('Bellows volt')
        ylabel('Flow [v/s]')
        %     volt_5pct=-prctile(-bellows_volt_drifted(:),5);
        %     volt_85pct=-prctile(-bellows_volt_drifted(:),85);
        
        %     if matlabpool('size') == 0 % checking to see if my pool is already open
        %         matlabpool 8
        %     end
        [~,~,grid_z]=ndgrid(ii,jj,kk);
        volt_ref_img=reshape(patient.bellows_volt_drifted(grid_z,patient.ref),size(ii,2),size(jj,2),size(kk,2));
        
        clear flow
        
        display('               Loading model parameters...');
        load([patient.model_folder '/model_params/alpha']);
        load([patient.model_folder '/model_params/beta']);
        load([patient.model_folder '/model_params/constant']);
        model_u=single([ constant.u(:) alpha.u(:) beta.u(:)]);% gamma_vec_u];
        model_v=single([ constant.v(:) alpha.v(:) beta.v(:)]);% gamma_vec_u];
        model_w= single([ constant.w(:) alpha.w(:) beta.w(:)]);% gamma_vec_u];
        clear alpha beta constant
        
        load([patient.folder '/hu_fit'])
        hu_fit_p1=squeeze(hu_fit(1,:,:,:));
        hu_fit_p2=squeeze(hu_fit(2,:,:,:));
        clear hu_fit
        
        for phase = 1:size(point,2)
            %     phase_pcntl=phases(phase)/100*(size(bv_smooth,1)+1);
            %     volt(phase)=bv_smooth(round(phase_pcntl))
            volt(phase) = breath_rep(point(phase));%-prctile(-bellows_volt_drifted(:),phases(phase));
            flow(phase) = flow_rep(point(phase));
            %     hu_correction=hu_fit_p2.*(volt_ref_img-volt);
            %     hu_correction=hu_fit_p2.*(volt(phase)-volt_ref_img);
            hu_fit_p2_finite=hu_fit_p2(isfinite(hu_fit_p2));
            hu_correction= nanmean(hu_fit_p2(patient.static_mask>0))*(volt(phase)-volt_ref_img);%ones(size(hu_fit_p2)).*
            
            figure(100);%plot([1,dim(3)],[volt(phase),volt(phase)],'r--','LineWidth',2);
            plot(volt(phase),  flow(phase),'yo','LineWidth',5)
            
            %         flow=0;
            
            
            E=[1 volt(phase) flow(phase)];
            
            vX = double(reshape(model_u*E', size(ii,2), size(jj,2), size(kk,2)));
            vY = double(reshape(model_v*E', size(ii,2), size(jj,2), size(kk,2)));
            vZ = double(reshape(model_w*E', size(ii,2), size(jj,2), size(kk,2)));
            
            %     dvf_pseudo = create_psuedo_dvf_vec_130215(folder,model_folder, static_mask, volt(phase),flow,scans,ii,jj,kk);
            
            %      cd(model_folder)
            %      mkdir('Psuedo_dvfs')
            %      cd('Psuedo_dvfs')
            %     save(sprintf('dvf_pseudo_phase_%d',phases(phase)),'dvf_pseudo')
            %     vX=squeeze(dvf_pseudo(:,:,1,:));
            %     vY=squeeze(dvf_pseudo(:,:,2,:));
            %     vZ=squeeze(dvf_pseudo(:,:,3,:));
            tic
            
        display(sprintf('Scan %d of %d; Inverting DVF (may take ~10mins)', phase,size(point,2)))
        [ivy,ivx,ivz]=backwards2forwards(vY,vX,vZ);
        ivy(ivy == 0) = NaN;
        ivx(ivx == 0) = NaN;
        ivz(ivz == 0) = NaN;
        
        display(sprintf('Deforming Scan %d of %d', phase,size(point,2)));
        image_ref = deform_image;%patient.static;%
         image_ref(patient.body_mask<1)= min(patient.static(:));
        image_scan_pseudo_tmp = deform_image_mirt3d(image_ref, ivy, ivx, ivz);
        image_scan_pseudo_tmp(isnan(image_scan_pseudo_tmp)) = min(image_scan_pseudo_tmp(:));
        mask_scan_pseudo = deform_image_mirt3d(patient.static_mask, ivy, ivx, ivz);
        mask_scan_pseudo(isnan(mask_scan_pseudo)) = 0;
%         mask_scan_pseudo(isnan(mask_scan_pseudo)) = min(mask_scan_pseudo(:));
        deformed_mask = logical(mask_scan_pseudo);
        %         deformed_image_tmp = deform_image_DT(img_average,vX, vY, vZ);
        %         deformed_mask_tmp= deform_image_DT(static_mask,vX, vY, vZ);
        %         deformed_mask=deformed_mask_tmp.image;
        %         deformed_mask(deformed_mask>0)=1;
        %         hu_correction(deformed_mask==0)=0;
        hu_correction(mask_scan_pseudo==0)=0;
        deformed_image=image_scan_pseudo_tmp+hu_correction;
%         [image_scan_pseudo, mask_scan_pseudo]  = deform_using_transformix_toolbox(patient, scan_folder);
        
%         image_scan = metaImageRead([patient.folder_static sprintf('/scan_%d_cut',phase)]);
            
            %
            %         display(sprintf('Deforming Image %d of %d', phase,num_phases));
            %         deformed_image_tmp = deform_image_DT(img_average,vX, vY, vZ);
            %         deformed_mask_tmp= deform_image_DT(static_mask,vX, vY, vZ);
            %         deformed_mask=deformed_mask_tmp.image;
            %         deformed_mask(deformed_mask>0)=1;
            %         hu_correction(deformed_mask==0)=0;
            %         deformed_hu_correction = deform_image_DT(hu_correction,vX, vY, vZ);
            %         hu_fit_p2_X0 = deform_image_DT(hu_fit_p2,vX, vY, vZ);
            %         hu_fit_p2_X0 = hu_fit_p2_X0.image;
            %         deformed_hu_correction = deformed_hu_correction.image;
            %         deformed_hu_correction(deformed_mask==0)=0
            %         deformed_image=deformed_image_tmp.image+deformed_hu_correction;
            %
            clear deformed_mask_tmp;
            clear deformed_image_tmp;
            
            %         partial_air=(1-((deformed_image+943)/(44+943)));
            %
            %         volume_simple(phase)=nansum(deformed_mask(:))*voxel_size*0.001/1000;
            %         volume_partial(phase)=nansum(partial_air(deformed_mask>0))*voxel_size*0.001/1000;
            
            save([patient.folder_4DCT  sprintf('/Deformed_pseudo_Image_phase_%d',phase)],'deformed_image')
            save([patient.folder_4DCT sprintf('/Deformed_pseudo_Mask_phase_%d',phase)],'deformed_mask')
            
            metaImageWrite(deformed_image,[patient.folder_4DCT  sprintf('/Deformed_pseudo_Image_phase_%d',phase)],'ElementSpacing', [1 1 1]);
            metaImageWrite(double(deformed_mask),[patient.folder_4DCT sprintf('/Deformed_pseudo_Mask_phase_%d',phase)],'ElementSpacing', [1 1 1]);

            %
            %      if phases(phase)==5;
            %          cd(model_folder)
            %          image_X0=deformed_image;
            %          mask_X0=deformed_mask;
            %          save('image_X0', 'image_X0')
            %          save('mask_X0', 'mask_X0')
            %          save('hu_fit_p2_X0', 'hu_fit_p2_X0')
            %      end
            toc
            
        end
        
        
    end
